
for j=flagstar-1:-1:1
%bs(j)=bs(j+1)-dbs(j+1)*ds;
s=ss(j);
betas=beta0+beta1*s;
lambdas=max(lambda0*(1-zeta0*s^zeta1),0);

zetas=max(eta0*(1-s^eta1),0);
zetan(j)=zetas;

ix=ixx(j);
x=xx(j);

for kkj=1:10000
db=(bs(j+1)-bs(j))/ds;

for kj=1:20
    


ixnew=(1-1/(bs(j)^(1-1/psi)*(A-ix-x)^(1/psi)/rho))/etaI;
ix=99/100*ix+ixnew/100;


end
errorode(j)=rho/(1-1/psi)*(((A-ix-x)/bs(j))^(1-1/psi)-1)+(ix-etaI/2*ix^2-deltaI)-gamma*sigma^2/2+lambdas/(1-gamma)*(betas/(betas+1-gamma)-1)+((x/s-etaX/2*x^2/s^2-deltaX)-(ix-etaI/2*ix^2-deltaI))*s*(bs(j+1)-bs(j))/ds/bs(j)+zetan(j)/(1-gamma)*((bsB(j)/bs(j))^(1-gamma)-1);
while abs(errorode(j))>10^(-12)
derrorode=rho/(1-1/psi)*(((A-ix-x)/bs(j))^(1-1/psi-1)*(1-1/psi)*(A-ix-x)/bs(j)^2*(-1))+((x/s-etaX/2*x^2/s^2-deltaX)-(ix-etaI/2*ix^2-deltaI))*s*((-1)/ds/bs(j)+(bs(j+1)-bs(j))/ds/bs(j)^2*(-1))+zetan(j)*(bsB(j)/bs(j))^(-gamma)*bsB(j)/bs(j)^2*(-1);
bs(j)=bs(j)-errorode(j)/derrorode;
errorode(j)=rho/(1-1/psi)*(((A-ix-x)/bs(j))^(1-1/psi)-1)+(ix-etaI/2*ix^2-deltaI)-gamma*sigma^2/2+lambdas/(1-gamma)*(betas/(betas+1-gamma)-1)+((x/s-etaX/2*x^2/s^2-deltaX)-(ix-etaI/2*ix^2-deltaI))*s*(bs(j+1)-bs(j))/ds/bs(j)+zetan(j)/(1-gamma)*((bsB(j)/bs(j))^(1-gamma)-1);
end
 xnew=(1-bs(j)/db*(1-etaI*ix))/etaX*s;
 x=999/1000*x+xnew/1000;
end

ixx(j)=ix;
xx(j)=x;

end






for j=1:flagstar-1
s=ss(j);
betas=beta0+beta1*s;
lambdas=max(lambda0*(1-zeta0*s^zeta1),0);
zetas=max(eta0*(1-s^eta1),0);
zetan(j)=zetas;

ix=ixx(j);
x=xx(j);    
    
db=(bs(j+1)-bs(j))/ds;
errori(j)=(1-1/(bs(j)^(1-1/psi)*(A-ixx(j)-xx(j))^(1/psi)/rho))/etaI-ixx(j);
%errorx(j)=(1-(bs(j)-s*db)/db*(1-etaI*ixx(j)))/etaX*ss(j)-xx(j);
errorx(j)=(1-bs(j)/db*(1-etaI*ixx(j)))/etaX*ss(j)-xx(j);
errorode(j)=rho/(1-1/psi)*(((A-ix-x)/bs(j))^(1-1/psi)-1)+(ix-etaI/2*ix^2-deltaI)-gamma*sigma^2/2+lambdas/(1-gamma)*(betas/(betas+1-gamma)-1)+((x/s-etaX/2*x^2/s^2-deltaX)-(ix-etaI/2*ix^2-deltaI))*s*(bs(j+1)-bs(j))/ds/bs(j)+zetan(j)/(1-gamma)*((bsB(j)/bs(j))^(1-gamma)-1);
end
errori(j+1)=0;
errorx(j+1)=0;
errorode(j+1)=0;

bsG=bs;


% 
% 
subplot(3,2,1)
plot(ss,bsG,'LineWidth',1.5)
hold on;
% 
% % plot(ss,bs0,'k:','LineWidth',1.5)
% % hold on;
% 
plot(ss,bsB,'r--','LineWidth',1.5)
hold on;




subplot(3,2,2)
plot(ss,ixx,'LineWidth',1.5)
hold on;

plot(ss,ixx0,'k:','LineWidth',1.5)
hold on;

subplot(3,2,3)
plot(ss,xx,'LineWidth',1.5)
hold on;

plot(ss,xx0,'k:','LineWidth',1.5)
hold on;


subplot(3,2,4)
plot(ss(1:flagstar),errorode(1:flagstar),'LineWidth',1.5)
hold on;
plot(ss(1:flagstar),errori(1:flagstar),'r--','LineWidth',1.5)
hold on;
plot(ss(1:flagstar),errorx(1:flagstar),'k:','LineWidth',1.5)
hold on;

subplot(3,2,5)
plot(ss,ixx-ixx0,'LineWidth',1.5)
hold on;

subplot(3,2,6)
plot(ss,bs-bs0,'LineWidth',1.5)
hold on;

% subplot(3,2,5)
% plot(ss,dbs-dbsan,'LineWidth',1.5)
% hold on;
